function [r2] = model_r2_mean_2D(t,L,mode,sigma,c2,Kr2_1,Kr2_2,Kr2_3,Kr2_4)
if nargin == 5
    [Kr2_1,Kr2_2,Kr2_3,Kr2_4] = lpcvd_Kr2_sincos(L,mode);
end

lambda = zeros(mode+1,mode+1);
for m = 0:mode
    for n = 0:mode
        lambda(m+1,n+1) =  -c2*4*(m^2+n^2);
    end
end

r2 = zeros(size(t));
z1 = zeros(mode+1,mode+1);
z2 = z1;
z3 = z1;
z4 = z1;
sigma2 = sigma^2;
for i = 1:length(t)
    for m = 0:mode
        for n = 0:mode
            if m == 0 && n== 0
                z2(m+1,n+1) = (sigma2)*t(i);
            elseif m ~= 0 && n~= 0
                z1(m+1,n+1) = sigma2*(exp(2*lambda(m+1,n+1)*t(i))-1)/(2*lambda(m+1,n+1));
                z2(m+1,n+1) = sigma2*(exp(2*lambda(m+1,n+1)*t(i))-1)/(2*lambda(m+1,n+1));
                z3(m+1,n+1) = sigma2*(exp(2*lambda(m+1,n+1)*t(i))-1)/(2*lambda(m+1,n+1));
                z4(m+1,n+1) = sigma2*(exp(2*lambda(m+1,n+1)*t(i))-1)/(2*lambda(m+1,n+1));
            elseif m == 0 && n~=0
                z2(m+1,n+1) = sigma2*(exp(2*lambda(m+1,n+1)*t(i))-1)/(2*lambda(m+1,n+1));
                z4(m+1,n+1) = sigma2*(exp(2*lambda(m+1,n+1)*t(i))-1)/(2*lambda(m+1,n+1));
            elseif m ~= 0 && n==0
%                 z1(m+1,n+1) = sigma2*(exp(2*lambda(m+1,n+1)*t(i))-1)/(2*lambda(m+1,n+1));
                z2(m+1,n+1) = sigma2*(exp(2*lambda(m+1,n+1)*t(i))-1)/(2*lambda(m+1,n+1));
                z3(m+1,n+1) = sigma2*(exp(2*lambda(m+1,n+1)*t(i))-1)/(2*lambda(m+1,n+1));
            else
                error('Error\n');
            end
        end
    end
    r2(i) = sum(sum(Kr2_1.*z1+Kr2_2.*z2+Kr2_3.*z3+Kr2_4.*z4));
%     r2(i) = sum(sum(Kr2_3.*z3));
end



